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Abstract. 

A system of atoms connected by harmonic springs to their nearest neighbors 
on a lattice is coupled to Ising spins that are in contact with a thermal bath and 
evolve under Glauber dynamics. Assuming a nearest-neighbor antiferromagnetic 
interaction between spins, we calculate analytically the equilibrium state. On a one- 
dimensional lattice, the system exhibits first and second order phase transitions. The 
order parameters are the total magnetization and the number of spin pairs in an 
antiferromagnetic configuration. On a hexagonal two dimensional lattice, spins interact 
with their nearest-neighbors and next-nearest-neighbors. Together with the coupling 
to atoms, these interactions produce a complex behavior that is displayed on a phase 
diagram. There are: ordered phases associated to ripples with atomic wavelength 
and antiferromagnetic order, ordered phases associated to ripples with nanometer 
wavelengths and ferromagnetic order, disordered glassy phases, and other phases 
presenting stripes formed by different domains. These static phases are discussed 
in relation to existing experiments and results for other models found in the literature. 
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1. Introduction 

Graphene is a one-atom thick crystal membrane with extraordinary mechanical and 
electronic properties pa ei ej. Experimental characterization of suspended graphene 
shows that it is covered with ripples. These ripples are several nanometers long waves 
of the sheet without a preferred direction PQE>], modify the electronic band structure [6], 
and are expected to be relevant in the understanding of electronic transport in graphene 
|7]. Also, more recently, buckling in which the unit structures consist of only two-three 
unit cells of the graphene honeycomb lattice has been experimentally observed [8j. 

There have been many studies of ripples. The earliest studies using Monte Carlo 
[9J or molecular dynamics simulations [_fO] have shown that ripples may be connected 
to variable length a bonds of carbon atoms and may be caused by thermal fluctuations. 
Other studies have explored the connection between rippling and electron-phonon 
coupling nn m- In particular, it has been suggested that, at zero temperature, 
the electron-phonon coupling may drive the graphene sheet into a quantum critical 
point characterized by the vanishing of the bending rigidity of the membrane [T3 j. The 
continuation of this work by J. Gonzalez [13] discusses the nonzero expectation value of 
the mean curvature (the Laplacian of the flexural phonon field) once the bending rigidity 
of the membrane vanishes, and its role as order parameter. Alternatively, assuming 
that the graphene membrane is fluctuating in 2 + d dimensions (with d > 1), Guinea 
et al. have calculated the dressed two-particle propagators of the elastic and electron 
interactions. They have found a collective mode which becomes unstable at a nonzero 
wave vector and causes the appearance of Gaussian curvature [15j. Amorim et al. [T6] 
estimate the crossover temperature between quantum and classical descriptions to be 
70-90 K. Thus a quantum description of ripples is not necessary at room temperature. 
All these studies investigate and characterize rippling as equilibrium phenomena. 

We are interested in rippling dynamics and stability of static corrugations under 
disturbances. In experiments to visualize ripples using a transmission electron 
microscope (TEM), the suspended graphene sheet is hit by a low-intensity electron 
beam that may push atoms out-of-plane upward or downward in a random fashion. 
An alternative technique to visualize ripples is using a scanning tunneling microscope 
(STM) [17], In this case, the graphene sheet is pushed and locally heated in the region 
close to the tip. Depending on the tunneling current and the voltage between tip and 
sheet, the latter may undergo a phase transition from a flexible (rippled) to a rigid 
(buckled) state [18] . 

In Ref. |19j . the authors simplify the distortion of the 2d crystal by modeling it with 
two-state spin-like variables (+1 upward, —1 downward). There are antiferromagnetic 
interactions among these spins, because the out-of-plane shift of the atoms in opposite 
directions stabilize the strictly 2d system while keeping the gapless band structure of 
graphene [19] - A rich phase diagram is found, including paramagnetic, ordered and 
glassy phases, depending on the temperature and the values of the nearest-neighbor 
and next-nearest-neighbor couplings. In this way, they describe the formation and 
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origin of the atomic scale rippling found in Ref. [8j. On the other hand, there are 
also models that investigate rippling by considering at each lattice site a continuous 
variable u describing the out-of-plane displacement of the carbon atom coupled to a 
spin variable (cr = ±1) representing an internal degree of freedom [20, 21]. This may be 
understood as a mechanical system coupled to a spin system. The spin at each lattice 
site represents the non-saturated fourth bond that, by a physical mechanism similar 
to the one discussed in [T9], tries to pull the corresponding atom upward (u > 0) or 
downward (u < 0) from the flat sheet configuration. Mathematically, this is done by 
introducing a linear coupling term in the system energy, proportional to —ucr for each 
lattice site. In these simple models, the mechanical system is either a chain of oscillators 
Bni or a discrete elasticity model of the hexagonal graphene lattice while the spins 
are in contact with a thermal bath and flip randomly according to Glauber dynamics at 
the temperature T of the thermal bath. In both models, the system forms metastable 
but long-lived ripples assuming slow spin relaxation [20, [2Tlj. 

ft is worth investigating a combination of the two approaches described in the 
previous paragraph. Firstly, it seems sensible to model the out-of-plane displacement 
at each lattice site by a continuous variable as in Refs. [20, El], which is driven by 
the internal degree of freedom represented by the spin. Secondly, these spins certainly 
interact among themselves, by the mechanism proposed in [T9] , Thus, in this paper 
we discuss the formation and dynamics of ripples in graphene through a system of 
atoms connected by harmonic springs and coupled to interacting Ising spins. We start 
from the spin-oscillator chain model [20| 2J. I22J and add interactions among spins that 
produce stable rippling states. There appear different phases and transitions between 
them depending on parameter values. 

The plan of the paper is as follows. In Section [2] we introduce the one-dimensional 
model and calculate analytically the equilibrium state. We also present the equations 
that determine the dynamics of the system. Section [3] is devoted to the expansion 
of the model to 2 dimensions on a hexagonal lattice, with first and second-neighbors 
interactions between spins. Besides, we do a systematic analysis of the system modifying 
the parameters and studying the stationary configuration after the transitory, using a 
phase diagram in Section 3.1, and describing the different phenomenology of each phase 
in Section 322 The main conclusions are presented in Sec. [4| Relevant information 
that is not covered in the main text is presented in the Appendices: some geometrical 
expressions for the hexagonal lattice are discussed in Appendix A while images of the 
different phases are collected in |Appendix B 


2. The one-dimensional model 

To start with, we consider a one-dimensional chain of N oscillators with nearest-neighbor 
interactions, in which each oscillator is linearly coupled to an Ising spin <jj = ±1. A 
detailed investigation of this model can be found in Ref. [20]. Therein, it was shown 
that, for appropriate temperatures, the stable equilibrium configuration has only one 
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ripple, although there appeared some more complex long-lived metastable rippled states. 
To explore whether stable multi-rippled equilibrium configurations are possible, we add 
an anti-ferromagnetic term to the hamiltonian, 


N 


« = £ 


P j k 


3=0 


2m ^ 2 f u 3 a j ^ <T j+3 <T j 


( 1 ) 


Here u and p are the vertical displacement and momentum respectively, and the extreme 
oscillators and spins are fixed (u 0 = p 0 = cr 0 = Un+i = Pn+i = &N+ 1 = 0) [23]. The 
dynamics of the model is as follows: (i) the oscillators’ equations of motion, 


milj — k (uj .|_i + Uj -1 — 2uj) = ftjj, 


( 2 ) 


are the usual ones, whereas the spins evolve according to Glauber dynamics [24]. The 
transition rate from the configuration (u,p,cr) to (u.p, Rjcr) (obtained from cr by 
flipping the j-th spin) is 


Pi 
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Wj(o-\ u ,p) = -(1 -PjCTj), 


tanh 


7 J, 

j u 3 ~ + a 3+ 1 ) 
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(4) 


in which a is the characteristic attempt rate for the spin flips. The Glauber transition 
rates ensure that the system satisfies detailed balance, and therefore the system 
reaches equilibrium for long enough times. In equilibrium, the probability of a certain 
configuration ( u,p,cr) is proportional to e~ n / T /Z, where we measure the temperature 
T in units of energy. 

As explained in |20j . for J = 0 rippling appears provided the temperature is less 


than 


T n = 


f 2 K. 


N 


K N 
K n ~ , 


, , -,v , (5) 

k 7r 

as N —» oo. To guarantee that the diffusive term in (J2]) remains finite in the continuum 
limit, it is convenient to nondimensionalize the equations of motion as follows: 


kun 


U] fK 2 N 


t* = 


t 


Kn V m’ 


J aKN\[rn 

K = —, 0 = - 7 =-, 

To Vk 


T 

6 = — — T 


k 


T n ' PKj ' 


L 0 J ^n 

Then the transition rates and the equations of motion become 
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We will omit the asterisks in nondimensional equations from now on. In order to 
obtain the scaling of the critical temperature, we need to know the scaling of the 
model parameters with the system size. In principle, this could be done by deriving 
our model from a fundamental microscopic one, but this is outside the scope of this 
paper. Nevertheless, we discuss some possible scalings in the following. If both the 
elastic constant k and the antiferromagnetic coupling constant J are considered to be 
independent of the system size, the only remaining parameter is /, the coupling between 
the elastic and internal (spin) degrees of freedom. If / is also independent of the system 
size, the critical temperature T 0 diverges as N 2 . In this case, rippling would be observed 
at all temperatures. On the other hand, a finite value of T 0 in the large system size 
limit is obtained when / oc N _1 . Then, rippling would be observed only for T < Tq. In 
principle, both situations are compatible with current experiments, in which rippling is 
observed over a wide temperature range. 

Let us consider now the equilibrium situation. Equation (|8]) can be averaged, with 
the result 

< 9 > 

in which ( u ) and (a) are the equilibrium average height and spin at position x = i/N\ 
the system has unit size in the continuous space variable x = i/N, 0 < x < 1. Therefore, 
the average curvature of the ripples is directly linked to the average spin. Very recently, 
this idea has been used to develop a phenomenological Ising model to study rippling in 
graphene in scanning tunneling microscopy experiments, in which each spin represents 
a whole ripple and the spin sign gives its corresponding convexity fl8]. Interestingly, we 
can derive an effective free energy for the string, by integrating the canonical distribution 
over the momenta p and the spins er: the resulting probability V becomes a functional 
of the string profile u(x), which in dimensionless variables reads [25] 



The quantity In f is the logarithm of the spins partition function per site, which depends 
on the “field” fuj/T —tu/9 and the coupling J/T —» k/6. The particularization of this 
free energy to the J = 0 case was obtained in Ref. M • For J jtz 0, in order to calculate 
In V the system is divided into a set of nearly independent subsystems with N s 1 
sites each, but such that N s N and the “field” u can be considered almost constant 
within each subsystem. We may denote the subsystems by S(x), in which x corresponds 
to the position of the subsystem in the continuum limit. In each subsystem, the local 
average magnetization p = N~ x YljeS(x)( a j) an d correlation C = N^ 1 Y^jeS(x)( a j a j+ 1 ) 
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as u plays the role of the external field and k is the coupling constant. Global order 
parameters may be defined as 


M = 




( 12 ) 


in which M is the absolute value of the total magnetization per site and the fraction 
of spin pairs in an antiferromagnetic configuration is VC = 1 (VC = 0) for perfect 
antiferromagnetic (ferromagnetic) ordering. 


Equation (10) clearly shows that the free energy is an extensive quantity if the 


dimensionless variables are of the order of unity, which is consistent with our scalings. 


Taking into account the expressions for the free energy, Eq. (10), and the average 


magnetization, Eq. the equation giving the equilibrium profile (J9]) is nothing but 
the Euler-Lagrange equation for the free energy functional, as it should be: the string 
profile minimizes the free energy. 

In what follows, we summarize the main physical implications of the short-ranged 
antiferromagnetic interaction, as compared to the J = 0 case. The flat profile u = 0 
is always a solution of Eq. & but it becomes unstable for exp(— 2 k/9)/9 > 1. For 
J = 0, there appear rippled configurations with non-zero magnetization that are stable 
for 9 < 1 BS|. For J ^ 0 the bifurcation condition exp(—2 k/9)/9 = 1 produces two 
temperatures 9\ and 9 2 for k < 0.18, as seen in the left panel of Fig. [TJ Specifically, 
these rippled configurations become unstable for low (high) enough temperatures, 6 < 6i 
(i 9 > 9 2 ). The transitions at 6\ and 9 2 are of second order, the order parameters M and 
VC are continuous because u bifurcates continuously from the solution u = 0, similarly 
to the behaviour found in the J = 0 case. For k > 0.18, this rippled ferromagnetic 
phase no longer exists because the antiferromagnetic coupling is too strong. 

In the limit 9 —> 0 + , the partition function of the spins becomes 


limdlnC Q, = k + (M - 2k)t7(|u| -2k), 


(13a) 


lim/r = sgn(u) rj(\u\ — 2k), lirnC' = sgn(|u| — 2k). (13b) 

where gfx) is the Heaviside step function, gfx) = 1 for x > 0 and gfx) = 0 for x < 0, 
and sgn(a:) is the sign function, sgn((c) = 2gfx) — 1. In the flat configuration, /i = 0 
and C — — 1 everywhere. There appears a new rippled low temperature phase, that 
is antiferromagnetic near the boundaries because of the clamped boundary conditions. 
Inside an interval of length x 0 close to the boundaries, |u < 2 k, and p, — 0, C — — 1 
therein. Then, the simplest configuration is composed of (i) two straight lines (u" = 0) 
near the boundaries, that is, in the intervals (0,x 0 ) and (1 — x 0 , 1), and (ii) a parabolic 
ripple with |w| > 2k ( u" = ±1) in between, for x G (x 0 ,1 — x 0 ), which corresponds to 
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Figure 1: (Left) Function controlling the bifurcation to the ferromagnetic phase. The 
bifurcation condition exp(— 2 k/6)/6 = 1 is equivalent to y{6 ) = — ^dlnd = k. It is 
clearly seen that the bifurcation condition is fulfilled by two temperatures 9 1 and 0 2 for 
k smaller than the maximum y ma . x = (2e) _1 of y , that is, k < y max . (Right) Function 
controlling the length Xq of the antiferromagnetic regions near the borders of the system 
at low temperatures. This length is given by 7t 2 xo( 1 — 2xo)/4 = k, which has two 
solutions Xqi and x 02 , %02 = 1/2 — x 0 i, for k < 7r 2 / 32 ~ 0.3 and no solutions for k > 0.3. 
At the limit value k = 0.3, it is Xqi = Xq 2 = 1/4. 


ferromagnetic ordering because y = ±1 and C — 1 for \u\ > 2 k. The continuity of u and 
u' at x = Xq implies that 7 t 2 o;o( 1 — 2x 0 ) = 4 k, which determines two possible values of x Q 
for k < 0.3, as seen in the right panel of Fig. [lj The configuration corresponding to the 
smallest value Xq\ is the absolute minimum of the free energy for 0 < Xq < 1/8, whereas 
both the flat string and the configuration corresponding to Xo 2 are metastable. For 
1/4 > xoi >1/8, the absolute minimum of the free energy corresponds to the flat string, 
whereas both configurations corresponding to Xoi an( 4 x 02 ar e metastable. The transition 
at x c = 1/8 (which corresponds to 25 percent of the spins being antiferromagnetic) is 
first order, because the order parameters M and T>C change discontinuously: in the zero 
temperature flat configuration, it is M = 0 and VC = 1, whereas in the rippled state we 
have M — 1 — 2xq and VC = 2xq. Moreover, as is usually the case in first-order phase 
transitions, the string has many other metastable configurations: they have n internal 
nodes Xi, i = 1,... ,n, at which u changes sign. The existence of the different phases 
and their relative stability will be thoroughly discussed elsewhere [25]. 

3. The two-dimensional model 

Here, we extend the model to dimension d— 2, in the hope that this will make it possible 
to find more complex behaviors. This extension is almost direct and, as we are interested 
in applying the model to mimic a graphene sheet, we consider a hexagonal lattice. Due 
to the symmetry, it is important to write the hamiltonian carefully. First, each atom is 
indexed: i will be the row index and j the column index, with the peculiarity that each 
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Figure 2: Figure summarizing the atoms indexes and the parameters of the unit cell of 
the hexagonal lattice. 


row comprises atoms with two different heights in a zigzag distribution, see Figure [2] 
It is important to note that the form of the equations will be qualitatively different for 
atoms for which \i — j\ is an even number (e-atoms), which have one nearest neighbor 
above and two below, and those for which \i — j\ is an odd number (o-atoms), which 
have one nearest neighbor below and two above, that is, the opposite situation. It is 
quite obvious that if the plane is rotated by an angle of 7r, the two types of atoms are 
interchanged. 

Taking into account the notation described above, we can write down the extension 
of the Id Hamiltonian to d — 2. Moreover, we introduce next-nearest-neighbor 
interactions, 


«=£ 


tj 


Pij / 

2^ f 'U'ijO'ij “t~ J “j - (Jij—2 "f" Ci- flj—l) 


+ 'y ^ j 2 \( U ij u i+l,j) + ( u ij u i,j-l) + ( u ij u i,j+l) ] 

\i— j’l^even 


dcqj + <7jj_i + &i : j + 1) 


(14) 


where i and j take values 1 —> i max and 1 —> j max , respectively. Following the same 
steps as in the previous section, the nondimensional equation of motion for each atom 
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(15) 

(16) 


7« = tanh[l^l - + 0ij+i + ffij-i) 

— ^( <7 ij'-2 + 0i,j+2 + 01-1,j-1 + 01-1,1+1 + 01+1,1-1 + 01+l,l+l)],(17) 

where A/v is a large scale parameter to be calculated later. As we said before, the 
difference between e-atoms with o-atoms follows from the rotation by 7r of the plane. 
For that reason, only equations for e-atoms have been written. In the latter equations, 
the height variable u and time are dimensionless. In the nondimensionalization, the 
same parameters as in equations (|6]) and ([7]) appear, with the addition of A = J'/T 0 , 
which corresponds to the new next-nearest-neighbor interaction. 

The length of each side of the finite hexagonal lattice is L = [3(n— 1) +1]Z/2, where 
l is the side of a unit hexagonal cell and n is the maximum value of the row index i in 
Fig. [ 2 J Let us measure all lengths in units of L, so that l — l/L tends to zero as the 
hexagonal lattice fills the plane. Then the expression within parenthesis in (15) has the 
limit |26l E3 Elj 


u i+ ij + Uij- 1 + u iJ+ 1 - 3 Uij — {d 2 x u + dyU), 

as a = \/3l —> 0. Therefore, we take Kjy proportional to a ', namely 


a/2 3n — 2 


Kn = —a 
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y/f 


OC 71, 


7T 


(18) 


(19) 


to guarantee that the diffusive term in (15) remains finite as l —> 0 (continuum limit). 
Note that the increments of the continuous variables are A x{i —> i + 1) = 31/2 and 
A y(j -+j + l) = a/2 = V^l/2, as seen in Fig. gl so that the hexagonal lattice goes to 


the unit square 0 < x, y < 1 in the continuum limit. For details, see Appendix A 


Once the system reaches the stationary state, equation (15) can be averaged 


ignoring thermal fluctuations. Thus, using equation (18) for n > 1 we get 

W), 


i v » = 


( 20 ) 


where ( u) and (a) are the average height and spin at the point (x, y) of the unit square. 
For k = A = 0, we have that (a) = tanh^/#) and the flat configuration ( u) = 0 
becomes unstable at 9 = 1, similarly to the situation in the Id case. This kind of rigidly 
buckled configurations have been observed in graphene in recent STM experiments [18j. 


Equation (20) tells us that there is a correspondence between lattice patterns given by 
the average height profile and the spin configuration. Specifically, the curvature of the 
rippling is directly proportional to the average spin. Therefore, in the following section 
we will mainly characterize the phases by the spin configuration. 
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3.1. Phase diagram 

As it has already been said, the stable steady state is a rigidly buckled configuration 
below the critical temperature (6 < 1), provided there is no interaction between spins. 
We expect that the introduction of the nearest neighbour interaction among the spins 
should introduce new phases. By analogy with the Id system, an antiferromagnetic 
nearest neighbour interaction should make antiferromagnetic ordered phases to appear 
for low enough temperatures. Looking for a more complex phenomenology, we introduce 


a next-nearest-neighbor interaction, as in Ref. |19j . This term appears in (14) through 
J' and in ( [T7| ) through its dimensionless counterpart A. 

ft is important to note that both J (k) and J' (A) may take positive or negative 
values, corresponding to antiferromagnetic and ferromagnetic interactions, respectively. 
However, only positive values of A will be considered, since a next-nearest-neighbor 
ferromagnetic interaction just strengthens the nearest-neighbor one [28]. For positive 
values of k and A the qualitative behavior is quite different. The nearest-neighbor 
interaction provides a defined minimum energy distribution in which each spin and 
its nearest-neighbors point in opposite directions. However, the next-nearest-neighbor 
interaction does not yield a defined minimum energy distribution. In fact, the second- 
neighbors of each atom are second-neighbors to each other, and therefore the next- 
nearest-neighbor interaction causes the system to be frustrated [H]. An enlightening 
discussion about frustration is given in the introduction of Eg. In principle, it is 
tempting to exclude negative values of k from the analysis. On intuitive grounds, one 
may conclude that the nearest-neighbour ferromagnetic coupling with k < 0 should only 
strengthen the already long-ranged ferromagnetic interaction among the spins induced 
by the spin-lattice coupling term —fuijOij m ■ Nevertheless, the situation is a little bit 
more complex, as discussed below. 

We plot a phase diagram to show in only one graph all the different behaviors, see 
Figure [3] In our simulations, we have chosen a nondimensional temperature 9 = 0.01, 
which is far below critical for k = A = 0. A key parameter is 


VC= — > 

N ^ 

\i—j\=even 


[3 + + CTij- 1 + Ojj'+l)], 


( 21 ) 


where N denotes the number of atoms in the lattice. This parameter estimates the 
domain-wall length |19j, and it is equal to 3 (resp. 0) for completely ferromagnetic 
(resp. antiferromagnetic) behavior. In addition, to delimit the regions on the diagram, 
we have used the absolute value of the usual magnetization 


M = 


v 


( 22 ) 


and energy fluctuations (proportional to the specific heat), 


F — y((A-e) 2 ). 


(23) 
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Figure 3: Phase diagram for a hexagonal lattice coupled to Ising spins. Different regions 
have been delimited using the domain-wall parameter, the magnetization and the specific 
heat. Once the equilibrium state is reached, each region has a different behavior, which 
is explained in the text. Also plotted is the line k/X — 4, which is a good estimate for 
the transition line between zones C and D. This agrees with the line separating phases 
Ordered 1 and 2 in Ref. [ T9] . 


Here e is the system energy and A*e = (e — (e)), where the angular brackets stands for 
the mean value that is calculated once the stationary state has been reached. 

3.2. Region characterization 

The different regions in the phase diagram have been characterized using the three 
parameters VC, M and F. Figure[3]is the superposition of the projections of M and VC 
on the plane X/6— tz/6. Each region of the plane correspond to different combinations of 
M and VC values. Once the regions have been delimited using the magnitudes described 
above, the system is allowed to evolve with k and A in one of the regions. Next, we 
verify that the system reaches the equilibrium state and we obtain the basic structures 
in the spin domains. Moreover, to check that we have actually reached the equilibrium 
state, another simulation is carried out with this distribution as the initial condition: if, 
aside from thermal fluctuations, no evolution is found, equilibrium has been reached. 

• Region A. VC ~ 3, M ~ 1. The plane is completely curved, and the spins are 
all pointing in the same direction. This situation corresponds to small values of k 
and A, for which the interaction that dominates is the one between the surface and 
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the spins, in agreement with the simple picture already present in the Id model, 
see Section [2] 

• Region B. It is the zone surrounding A, on which VC and M decrease from the 
A values to those on the other regions. Here, the system displays a behavior that 
is analogous to the one described at the end of Section [2] (in Id). The interaction 
between the surface and the spins is an effective ferromagnetic interaction with an 
intensity that decreases from the center to the border. Thus, the plane is curved 
but the spins close to the border are antifcrromagnetically arranged. 

• Region C. VC ~ 0.5, M ~ 0. The predominant interaction is the antiferromagnetic 
first-neighbor one. The equilibrium state (starting from a random initial spin 
distribution) is composed of antiferromagnetic domains. 

• Region D. VC increases from 0.5 to 1.2, M ~ 0. The states in this region are 
metastable. Taking the distribution corresponding to the equilibrium state of C or 
E as the initial condition, the system does not evolve to the other state, at least 
in a simulation time much greater than the relaxation time from random initial 
conditions. 

• Region E. VC ~ 1.2, M ~ 0. In this case, the spins are distributed in rows 
of two atoms in the lowest energy configuration. Beginning with random initial 
conditions, these two-atoms domains were created, with the rows in any of the 
three symmetrical directions. 

• Region F. VC ~ 1.5, M ~ 0.2. The interaction between the plane and the spins is 
relevant again, since the antiferromagnetic next-nearest-neighbor interaction (with 
k ~ 0) has no defined minimum energy distribution. The plane is curved, leading 
to a non-zero magnetization. 

• Region G. VC ~ 1.8, M ~ 0. The typical equilibrium configurations are long 
serpentine lines, with zero magnetization. Taking as initial conditions the spins 
arranged in rows, the system remains static. 

• Regions H and /. In them, the system evolves from the characteristic configurations 
of G to the ferromagnetic configurations of J, with the difference that in E[ the 
magnetization is different from zero whereas in / it is not. / is a ferromagnetic 
first-neighbor state, but with domains smaller than in J (VC smaller than in J). 

• Region J. VC ~ 2.5, M ~ 0.1. In this case the system behaves as a completely 
ferromagnetic first-neighbor system. Starting from random initial conditions, 
ferromagnetic domains grow until reaching a stationary state. In this case M is 
close to zero since spins are pointing to different directions in adjoining domains. 


The plots of the typical equilibrium configurations for each region are in Appendix B| 
It should be noted that our phase diagram does not contain a paramagnetic state 
because the chosen temperature, 6 = 0.01, is far below the critical temperature for k = 0 
and A = 0 (unity in our dimensionless variables). Each point of the phase diagram 
corresponds to the energy minimum to which the system evolves for the considered 
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parameters. Once it is in the neighborhood of this minimum, the energy barriers are so 
high that ergodicity is no longer valid, and the system remains frozen [30]. This causes 
an Edwards-Anderson order parameter 135 , 


Qea 


1 

N 


ij 


Rij — ij)i 


(24) 


to be different from zero at every point of the plotted phase diagram. On the other hand, 
close to k = X = 0, the order parameter qEA will vanish as 0 -A- oo, once ergodicity 
is recovered. In Eq. (24) the average should be understood as a time average or an 
extended Gibbs average in a phase space composed of disjoint ergodic components ffl- 


4. Conclusions 


We have studied a system of atoms connected by harmonic springs and coupled to 
Glauber spins. The spins are in contact with a thermal bath and interact with their 
neighbors. The Id system forms one ripple and becomes antiferromagnetic at the 
boundaries as p increases, until it becomes completely antiferromagnetic. When the 
system is on a 2d hexagonal lattice, each spin interacts with its nearest-neighbors and 
next-nearest-neighbors, aside from the coupling with the out-of-plane displacement. 
This situation generates different phases which are included in a phase diagram. 

The range of parameters in our phase diagram includes negative values for the 
nearest neighbor coupling constant k and is thus wider than the one used in [T5], in 
which only antiferromagnetic interactions were considered. The change in the sign of 
the spin-spin interaction can be produced by the scattering of the conduction electrons 
at the spins, see [321 133. We are interested in zero magnetization phases since they 
correspond to no overall bending. Our model provides different phases obeying this 
constraint: I and J are long wave length phases, similar to those observed in M, 
whereas C, E and G are phases with atomic wave length. G is a stripy phase (see Figure 


B1), which could be associated with patterns seen in |3]. The atomic wave length phases 
C and E correspond with phases Ordered 1 and 2, respectively, from Ref. ra- Therein, 
the line between these two phases is (in our variables) k/X = 4, which agrees with the 
limit of true stability of C here. Interestingly, neither the metastablc phase D or the 
other phases (including the long wavelength phases I and J) were found in Ref. [19]. 
In that reference, (i) only positive values of n were considered, and (ii) there was no 
spin-atom coupling. 

The buckling phase A is surrounded by rippled phases C, E, and G with no overall 
bending. Starting from a point of the phase diagram belonging to region C (rippled 
phase), if we increase the temperature while keeping k and A (supposed temperature 
independent) fixed, we move along a straight line of slope X/n in Figure [3] from phase C 
to the buckled phase A. In experiments with STM at fixed current, the temperature is 
locally increased at the tip region and this triggers a transition from a rippled flexible 
phase to a rigid buckled phase [18]. Thus our model contains the ripples-to-buckling 
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transition observed in experiments although more work needs to be done to explain 
STM observations in detail [25] . 

To conclude, our model is based in a few parameters controlling simple interactions 
which generate complex collective behaviors. This allows us to identify the interactions 
responsible for each pattern. In addition, the elastic feature of the model makes 
it possible to visualize and quantify the magnitude of the rippling, which could be 
compared with experiments once height measurements had been improved. 
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Appendix A. Geometrical expressions 


We want our hexagonal lattice to have equal overall length and height. Let n = i max be 
the total number of rows. Then 


] m ax = IntegerPart 


3 (n - 1) + 1 

71 


+ 1 


(A.l) 


is the total number of columns, and the height of the hexagonal lattice is 

3(n — 1) + P 


L = 


-l, 


(A.2) 


where l is the length of the side of a unit hexagonal cell. With these expressions, if 
n = 25, then j m ax — 43 and the vertical and horizontal side of the lattice are 1 and 
0.997 respectively, in units of L. Our finite hexagonal lattice is then roughly inscribed 
in a square and the nondimensional side of the unit hexagonal cell is l = l/L = 0.027. 


Appendix B. Phase diagram images 

In our simulations, we have used a lattice of 2,150 atoms and a temperature 6 = 0.01. 
We need to impose initial and boundary conditions for the membrane and the spins. 
Initially, the spins are in a completely random configuration, whereas the membrane is 
flat and at rest. The membrane is clamped (zero displacement) at the boundaries. As 
nearest and next-nearest neighbors determine the dynamics of a given spin, see eq. 
a spin located next to the boundary condition needs data from nearby spins located at 
the boundaries and also outside the lattice. The simplest possibility is that the spins of 
clamped boundary atoms and their nearest neighbors outside the lattice do not interact 
with the others, which can be achieved by formally assigning spin zero to them. 
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Figure Bl: Final configuration of the plane and the spins (red for spin up and blue for 
spin down) for different values of k and A, corresponding to different regions of the phase 
diagram in Fig. [3j From top to bottom and left to right: Region A, k/6 = —10 and 
X/9 = 2, Region B , k/6 = 0 and X/6 = 10, Region C, k/6 = 40 and X/6 = 3, Region E, 
k/6 = 40 and X/6 = 27, Region F, k/6 = 0 and X/6 = 24, Region G, k/6 = —25 and 
X/6 = 24, Region H, k/6 = —60 and X/6 = 30, Region /, k/6 = —73 and X/6 = 27, 
Region J, k/6 = —91 and X/6 = 3. 
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